#include "Mesh.h"

#include <cfloat>
#include <cassert>
#include <cmath>
#include <algorithm>

#define ACCURACY 1e-7

using namespace Mesh;

// 基本运算
vector3f vector3f::operator*(double t)
{
	return {x * t, y * t, z * t};
}

vector3f vector3f::operator/(double t)
{
	assert(fabs(t) > ACCURACY);
	return {x / t, y / t, z / t};
}

vector3f vector3f::operator-(const vector3f &v)
{
	return {x - v.x, y - v.y, z - v.z};
}

vector3f vector3f::operator+(const vector3f &v)
{
	return {x + v.x, y + v.y, z + v.z};
}

vector3f Mesh::operator-(vector3f &v)
{
	return {-v.x, -v.y, -v.z};
}

vector3f point3f::operator-(const point3f &p)
{
	return {x - p.x, y - p.y, z - p.z};
}

// 基本运算
coord4f coord4f::operator*(double t)
{
	coord4f c = *this;
	if (w == 1)
	{
		c.x /= t;
		c.y /= t;
		c.z /= t;
	}
	else
	{
		c.x *= t;
		c.y *= t;
		c.z *= t;
	}
	return c;
}

coord4f coord4f::operator/(double t)
{
	assert(fabs(t) > ACCURACY);
	coord4f c = *this;
	if (w == 1)
	{
		c.x *= t;
		c.y *= t;
		c.z *= t;
	}
	else
	{
		c.x /= t;
		c.y /= t;
		c.z /= t;
	}
	return c;
}

coord4f coord4f::operator-(const coord4f &v)
{
	return {x - v.x, y - v.y, z - v.z, w - v.w};
}

coord4f coord4f::operator+(const coord4f &v)
{
	coord4f c = {x + v.x, y + v.y, z + v.z, w + v.w};
	if (c.w == 2)
	{
		c.x /= c.w;
		c.y /= c.w;
		c.z /= c.w;
	}
	return c;
}

coord4f Mesh::operator-(coord4f &v)
{
	return {-v.x, -v.y, -v.z, v.w};
}

// 颜色直接叠加
rgb_space rgb_space::operator+(const rgb_space &c1)
{
	rgb_space c = {c1.r + r, c1.g + g, c1.b + b};
	c.r = (c.r > 1) ? 1 : c.r;
	c.g = (c.g > 1) ? 1 : c.g;
	c.b = (c.b > 1) ? 1 : c.b;
	return c;
}

// 最后释放内存
HEMesh::~HEMesh()
{
	for (auto it = m_vList.begin(); it != m_vList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
	for (auto it = m_eList.begin(); it != m_eList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
	for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
		if ((*it) != nullptr)
			delete *it;
}

// 创建顶点，并返回指针
vertex3f *HEMesh::CreateVertex(point3f p)
{
	vertex3f *vp = new vertex3f({p, nullptr});
	m_vList.push_back(vp);
	return vp;
}

// 创建三角面，并自动衔接
face *HEMesh::CreateFace(vertex3f *v0, vertex3f *v1, vertex3f *v2)
{
	// 判断三点是否共线
	if (onLine(v0->pos, v1->pos, v2->pos))
		return nullptr;

	face *f = new face;

	// 创建三条半边
	half_edge *e = new half_edge;
	half_edge *ne = new half_edge;
	half_edge *pe = new half_edge;
	m_eList.push_back(e);
	m_eList.push_back(ne);
	m_eList.push_back(pe);

	// 面持有任意一个半边
	f->e = e;

	// 初始化半边
	e->v = v0;
	e->incFace = f;
	e->next = ne;
	e->prev = pe;
	e->twin = nullptr;

	ne->v = v1;
	ne->incFace = f;
	ne->next = pe;
	ne->prev = e;
	ne->twin = nullptr;

	pe->v = v2;
	pe->incFace = f;
	pe->next = e;
	pe->prev = ne;
	pe->twin = nullptr;

	// 如果顶点没有半边，就赋予半边
	if (v0->e == nullptr)
		v0->e = e;
	if (v1->e == nullptr)
		v1->e = ne;
	if (v2->e == nullptr)
		v2->e = pe;

	// 遍历整个面表，逐一尝试建立 twin 连接
	// 优化：不再遍历面表，而是选择与顶点相邻的面尝试建立 twin 连接；当面的数量较大时具有优势
	// 再优化：由于三角形顶点两两共用面，只需要计算两个顶点的邻面
	// 注意要先取出所有邻面再建立连接，防止出现与自己连接的情况
	std::deque<face *> fList0 = adjacentFace(v0);
	std::deque<face *> fList1 = adjacentFace(v1);
	for (auto it = fList0.begin(); it != fList0.end(); ++it)
		stitchFace(f, *it);
	for (auto it = fList1.begin(); it != fList1.end(); ++it)
		stitchFace(f, *it);

	// 最后推入面
	m_fList.push_back(f);

	return f;
}

// 删除顶点
void HEMesh::DeleteVertex(vertex3f *v)
{
	// 获得所有相邻的面，将它们删除
	std::deque<face *> fList = adjacentFace(v);
	for (auto it = fList.begin(); it != fList.end(); ++it)
		DeleteFace(*it);

	// 最后删除顶点
	m_vList.erase(std::find(m_vList.begin(), m_vList.end(), v));
	delete v;
	v = nullptr;
}

// 删除面
void HEMesh::DeleteFace(face *f)
{
	// 获得三条半边
	half_edge *e[] = {f->e, f->e->next, f->e->prev};

	// 调整每条半边起点保存的半边
	for (int i = 0; i < 3; i++)
	{
		// 获得半边的起源点，和以 v 为起点的半边
		vertex3f *v = e[i]->v;

		// 删除面时，如果面上的顶点孤立，说明出现了非流形网格
		// 这时候其实不用考虑这个顶点了
		std::deque<half_edge *> eList = adjacentHalfEdge(v);

		// 如果 v 不止一条半边，就赋予新的半边，否则 v 将会孤立
		for (auto it = eList.begin(); it != eList.end(); ++it)
		{
			// 赋予不同的半边
			if ((*it) != e[i])
			{
				v->e = *it;
				break;
			}
		}
	}

	// 删除每条半边
	for (int i = 0; i < 3; i++)
	{
		// 首先将每条半边的反向边的反向边设为空
		half_edge *te = e[i]->twin;
		if (te != nullptr)
			te->twin = nullptr;

		// 如果此半边的起点仍然持有此半边，那么起点将会孤立，将该点保有的半边设为空
		if (e[i]->v->e == e[i])
		{
			e[i]->v->e = nullptr;
			e[i]->v = nullptr;
		}

		// 然后就可以删除此半边
		m_eList.erase(std::find(m_eList.begin(), m_eList.end(), e[i]));
		delete e[i];
		e[i] = nullptr;
	}

	// 最后删除面
	m_fList.erase(std::find(m_fList.begin(), m_fList.end(), f));
	delete f;
	f = nullptr;
}

// 清除孤立点
void HEMesh::DeleteDiscrete()
{
	// 找到所有孤立点
	for (auto it = m_vList.begin(); it != m_vList.end();)
	{
		// 如果不存在以它为起点的半边，说明此顶点孤立，直接删除
		if ((*it)->e == nullptr)
			it = m_vList.erase(it);
		else
			++it;
	}
}

// 是否有法向
bool HEMesh::HasVertexNormal() const
{
	return !m_nList.empty();
}

// 是否有颜色
bool HEMesh::HasVertexColor() const
{
	return !m_cList.empty();
}

// 判断网格是否是凸包
bool HEMesh::IsConvexHull() const
{
	// 对每个顶点，遍历每一个面
	for (auto vit = m_vList.begin(); vit != m_vList.end(); ++vit)
	{
		// 要排除离散点
		if ((*vit)->e == nullptr)
			continue;

		// 判断顶点是不是在所有面的上方，如果在某个面的上面，说明不是凸包
		for (auto fit = m_fList.begin(); fit != m_fList.end(); ++fit)
			if (aboveFace((*vit)->pos, *fit))
				return false;
	}

	return true;
}

// 判断点是否在网格的凸包中
bool HEMesh::IsInMesh(point3f p) const
{
	// 遍历每一个面，p 如果在某个面的上面，说明 p 在网格外
	for (auto it = m_fList.begin(); it != m_fList.end(); ++it)
		if (aboveFace(p, *it))
			return false;
	return true;
}

// 颜色线性插值
rgb_space Mesh::linear(rgb_space c1, rgb_space c2, double t)
{
	return {c1.r * (1 - t) + c2.r * t,
			c1.g * (1 - t) + c2.g * t,
			c1.b * (1 - t) + c2.b * t};
}

// 点线性插值
point3f Mesh::linear(point3f p1, point3f p2, double t, double a, double b)
{
	return {(p1.x * (b - t) + p2.x * (t - a)) / (b - a),
			(p1.y * (b - t) + p2.y * (t - a)) / (b - a),
			(p1.z * (b - t) + p2.z * (t - a)) / (b - a)};
}

// 计算点的距离
double Mesh::distance(point3f p1, point3f p2)
{
	return sqrt(pow(p1.x - p2.x, 2) + pow(p1.y - p2.y, 2) + pow(p1.z - p2.z, 2));
}

// 行列式 (v1 x v2) * v3
double Mesh::det(vector3f v1, vector3f v2, vector3f v3)
{
	return dot(cross(v1, v2), v3);
}

// 向量点积
double Mesh::dot(vector3f v1, vector3f v2)
{
	return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}

// 向量叉积
vector3f Mesh::cross(vector3f v1, vector3f v2)
{
	return {v1.y * v2.z - v1.z * v2.y,
			-v1.x * v2.z + v1.z * v2.x,
			v1.x * v2.y - v1.y * v2.x};
}

// 向量规范化
vector3f Mesh::normalize(vector3f p)
{
	double norm = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
	p.x /= norm;
	p.y /= norm;
	p.z /= norm;
	return p;
}

// 向量长度
double Mesh::length(vector3f p)
{
	return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
}

// 判断三点是否共线
bool Mesh::onLine(point3f p1, point3f p2, point3f p3)
{
	vector3f vec1 = p1 - p2;
	vector3f vec2 = p1 - p3;

	// 叉积足够大，则认为不共线
	if (length(cross(vec1, vec2)) > ACCURACY)
		return false;

	return true;
}

// “缝合”两个面，返回是否成功
bool Mesh::stitchFace(face *f1, face *f2)
{
	// 先取得公共半边
	std::deque<half_edge *> e = commonOnFace(f1, f2);

	// 如果有公共边就连接
	if (e[0] != nullptr && e[1] != nullptr)
	{
		// 如果公共边已经建立，但是不是两个面的公共边，则返回 false
		if ((e[0]->twin == nullptr || e[0]->twin->twin == e[0]) &&
			(e[1]->twin == nullptr || e[1]->twin->twin == e[1]))
		{
			// 然后互相为 twin
			e[0]->twin = e[1];
			e[1]->twin = e[0];
			return true;
		}
	}
	// 报错说明这个面破坏了网格的流形结构，或者没有公共边
	return false;
}

// 判断一个点是否在面上
bool Mesh::onFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 找到面的一个顶点
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 如果与法向正交，说明共面
	return (fabs(dot(v, faceNormal(f))) < ACCURACY) ? true : false;
}

// 判断一个点是否在面的上方
bool Mesh::aboveFace(point3f p, face *f)
{
	// 排除面是空的情况
	if (f == nullptr)
		return false;

	// 取一个顶点计算夹角
	point3f p0 = f->e->v->pos;
	vector3f v = p - p0;

	// 这里精度需要仔细考虑，不然判断会有问题
	return (dot(v, faceNormal(f)) > ACCURACY) ? true : false;
}

// Gouraud 插值：计算顶点的平均法向
vector3f Mesh::Gouraud(vertex3f *v)
{
	// 先获取所有相邻的面
	std::deque<face *> fList = adjacentFace(v);

	// 累加然后规范化
	vector3f n = {0, 0, 0};
	for (auto it = fList.begin(); it != fList.end(); ++it)
		n = n + faceNormal(*it);
	return normalize(n);
}

// 计算单位面法向
vector3f Mesh::faceNormal(face *f)
{
	// 找到面的三个顶点计算法向
	point3f v0 = f->e->v->pos;
	point3f v1 = f->e->next->v->pos;
	point3f v2 = f->e->next->next->v->pos;
	return normalize(cross(v1 - v0, v2 - v1));
}

// 获得边的目标点
vertex3f *Mesh::destination(half_edge *e)
{
	if (e == nullptr)
		return nullptr;
	e = e->next;
	return e->v;
}

// 获得一个面上的所有顶点
std::deque<vertex3f *> Mesh::vertexOnFace(face *f)
{
	std::deque<vertex3f *> vList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		vList.push_back(e->v);
		e = e->next;
	}
	return vList;
}

// 获得一个面上的所有半边
std::deque<half_edge *> Mesh::halfEdgeOnFace(face *f)
{
	std::deque<half_edge *> eList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		eList.push_back(e);
		e = e->next;
	}
	return eList;
}

// 获取两个面的公共半边
std::deque<half_edge *> Mesh::commonOnFace(face *f0, face *f1)
{
	std::deque<half_edge *> eList(2, nullptr);

	half_edge *e0 = f0->e;
	half_edge *e1 = f1->e;

	// 遍历两个面上的所有边
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			// e0 的起点是 e1 的终点，e1 的起点是 e0 的终点时，说明是公共边
			// 储存后直接返回
			if (e0->v == destination(e1) && e1->v == destination(e0))
			{
				eList[0] = e0;
				eList[1] = e1;
				return eList;
			}
			e1 = e1->next;
		}
		// 注意 e1 返回原位
		e1 = f1->e;
		e0 = e0->next;
	}
	return eList;
}

// 获得相邻的所有顶点
std::deque<vertex3f *> Mesh::adjacentVertex(vertex3f *v)
{
	std::deque<vertex3f *> vList;

	// 先获得以 v 为起点的半边，其目标点就是相邻顶点
	std::deque<half_edge *> eList = adjacentHalfEdge(v);
	for (int i = 0; i < eList.size(); i++)
		vList.push_back(destination(eList[i]));

	return vList;
}

// 获得以 v 为起点的半边
std::deque<half_edge *> Mesh::adjacentHalfEdge(vertex3f *v)
{
	std::deque<half_edge *> eList;
	half_edge *e = v->e;

	// 先顺时针寻找
	while (e != nullptr)
	{
		eList.push_back(e);
		e = e->twin;

		// 反向边，然后找到下一条边
		if (e != nullptr)
			e = e->next;

		// 如果回到初始边，就退出
		if (e == v->e)
			break;
	}

	// 如果没有回到初始边，就逆时针寻找，排除孤立顶点
	if (e == nullptr && v->e != nullptr)
	{
		e = v->e->prev->twin;
		while (e != nullptr)
		{
			eList.push_back(e);
			e = e->prev->twin;
		}
		// 这里不需要判断回到初始边，因为如果一个方向不能返回，那么反过来也不能返回
	}
	return eList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::adjacentFace(vertex3f *v)
{
	// 先获得以 v 为起点的半边
	std::deque<half_edge *> eList = adjacentHalfEdge(v);

	// 每条半边左边的面就是相邻的面
	std::deque<face *> fList;
	for (auto it = eList.begin(); it != eList.end(); ++it)
		fList.push_back((*it)->incFace);

	return fList;
}

// 获得相邻的所有面
std::deque<face *> Mesh::adjacentFace(face *f)
{
	std::deque<face *> fList;
	half_edge *e = f->e;
	for (int i = 0; i < 3; i++)
	{
		if (e->twin != nullptr)
			fList.push_back(e->twin->incFace);
		e = e->next;
	}
	return fList;
}

// 尝试用网格点创建四面体，返回是否成功
bool Mesh::tetrahedron(HEMesh *mesh)
{
	// 总顶点数至少是 4，并且还没有创建任何面和边
	if (mesh->m_vList.size() < 4 || mesh->m_fList.size() > 0 || mesh->m_eList.size() > 0)
		return false;

	// 先用前 3 个顶点建立一个三角形
	auto it = mesh->m_vList.begin();
	vertex3f *v1 = *it++;
	vertex3f *v2 = *it++;
	vertex3f *v3 = nullptr;
	vertex3f *v4 = nullptr;

	// 第三个顶点不能与前两个共线
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共线的顶点，方便起见，将它插入到头部
		if (!onLine(v1->pos, v2->pos, (*it)->pos))
		{
			v3 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v3);
			break;
		}
	}
	// 如果没找到，说明全共线，退出
	if (v3 == nullptr)
		return false;

	// 创建三角形
	face *f = mesh->CreateFace(v1, v2, v3);

	// 寻找第 4 个顶点，不与前三个共面
	it = mesh->m_vList.begin();
	it++, it++, it++;
	for (; it != mesh->m_vList.end(); ++it)
	{
		// 找到不共面的顶点，方便起见，将它插入到头部
		if (!onFace((*it)->pos, f))
		{
			v4 = *it;
			mesh->m_vList.erase(it);
			mesh->m_vList.push_front(v4);
			break;
		}
	}
	// 如果没找到，说明全共面，退出
	if (v4 == nullptr)
		return false;

	// 检查 4 顶点是否在面的上方
	bool above = aboveFace(v4->pos, f);
	if (above)
	{
		// 如果在上方，就要删除原来的面，将顺序反转，以确保法向正确
		mesh->DeleteFace(f);
		mesh->CreateFace(v3, v2, v1);
		mesh->CreateFace(v1, v4, v3);
		mesh->CreateFace(v3, v4, v2);
		mesh->CreateFace(v1, v2, v4);
	}
	else
	{
		mesh->CreateFace(v1, v3, v4);
		mesh->CreateFace(v3, v2, v4);
		mesh->CreateFace(v1, v4, v2);
	}

	return true;
}

// 取网格点的凸包
bool Mesh::convexHull(HEMesh *mesh)
{
	// 先尝试创建一个四面体
	if (!tetrahedron(mesh))
		return false;

	// 增量算法，对每个四面体以外的顶点，计算被这个顶点照射到的面
	// 将这些面删除，然后把面的边与顶点建立新的面
	auto vit = mesh->m_vList.begin();
	vit++, vit++, vit++, vit++;
	for (; vit != mesh->m_vList.end(); ++vit)
	{
		// 标记被照射的面，只需要保存一个面的指针，因为所有标记面都相邻
		std::deque<face *> fList;
		for (auto fit = mesh->m_fList.begin(); fit != mesh->m_fList.end(); ++fit)
		{
			// 在表面上面，说明被照射到，将其记录下来
			if (aboveFace((*vit)->pos, *fit))
			{
				(*fit)->identity = 1;
				fList.push_back(*fit);
			}
			else
				(*fit)->identity = 0;
		}

		// 检查被照射的面，标记相邻的未照射面
		std::deque<half_edge *> eList;
		for (auto fit = fList.begin(); fit != fList.end(); fit++)
		{
			// 这里不需要考虑没有反向边，因为凸包的边一定有反向边
			half_edge *e = (*fit)->e;
			half_edge *ne = e->next;
			half_edge *pe = e->prev;

			// 记录所有反向边
			if (e->twin->incFace->identity != 1)
				eList.push_back(e->twin);
			if (ne->twin->incFace->identity != 1)
				eList.push_back(ne->twin);
			if (pe->twin->incFace->identity != 1)
				eList.push_back(pe->twin);
		}

		// 删除光照面。直接删除不会有问题，如果出现顶点指向的边为空的情况
		// 有两种可能：1. 离散点，则在最后会被删除；2. 处于被光照和未被光照的边界上，则新创建的面会重新给该顶点赋予一个半边
		for (auto it = fList.begin(); it != fList.end(); ++it)
			mesh->DeleteFace(*it);

		// 遍历所有反向边，按照连接顺序，每次要创建具有公共边的新面，防止出现非流形网格
		// 这样创建的面都会相邻，不会出现恰好有一个公共顶点的两个面，而周围都是空的情况
		int count = eList.size();
		auto eit = eList.begin();
		while (count > 0)
		{
			mesh->CreateFace(destination((*eit)), (*eit)->v, *vit);
			count--;

			// 找到目标点对应的边
			for (auto it = eList.begin(); it != eList.end(); ++it)
			{
				if ((*it)->v == destination((*eit)))
				{
					eit = it;
					break;
				}
			}
		}
	}

	// 最后移除所有离散点
	mesh->DeleteDiscrete();

	return true;
}

// 修复网格，有时删除操作造成某些非离散点存放的半边为空，此函数修复这一问题，返回是否进行了修复
bool Mesh::fixMesh(HEMesh *mesh)
{
	bool isBad = false;
	for (auto it = mesh->m_fList.begin(); it != mesh->m_fList.end(); ++it)
	{
		// 获得三条半边
		half_edge *e[] = {(*it)->e, (*it)->e->next, (*it)->e->prev};
		for (int i = 0; i < 3; i++)
		{
			if (e[i]->v->e == nullptr)
			{
				e[i]->v->e = e[i];
				isBad = true;
			}
		}
	}
	return isBad;
}